#include <iostream>
#include <iomanip>
#include <fstream>

#include "TFile.h"
#include "TString.h"
#include "TTree.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "TPad.h"
#include "TGraph.h"
#include "TStyle.h"

#include "CommandLineInterface.hh"
#include "Peaks.hh"
using namespace std;


const char* detector[] = {"Forward_gam_pad_0","Forward_gam_pad_1","Forward_gam_pad_2","Forward_gam_pad_3","Backward_gam_pad_0","Backward_gam_pad_1","Backward_gam_pad_2","Backward_gam_pad_3"};


int lookup(char det[]){
  const int nrdet = sizeof(detector)/sizeof(*detector);
 
  for(int k=0; k!= nrdet; k++)
    if( strcmp(det,detector[k]) == 0 )
      return k;
  return -1;
}

int main(int argc, char* argv[]){
  char* InputFile = NULL;
  char* OutputFile = NULL;
  char* SettingFile = NULL;
  char* HOutputFile = NULL;
  char* Name = NULL;
  CommandLineInterface* interface = new CommandLineInterface();

  interface->Add("-i", "inputfile", &InputFile);
  interface->Add("-o", "outputfile", &OutputFile);
  interface->Add("-s", "settingsfile", &SettingFile);
  interface->Add("-h", "also produce histogramm output file _hist.root", &HOutputFile);
  interface->CheckFlags(argc, argv);

  if(InputFile == NULL || OutputFile == NULL || SettingFile == NULL ){
    cerr<<"You have to provide input, output and settings file!"<<endl;
    exit(1);
  }
  cout<<"input file: "<<InputFile<< endl;
  cout<<"output file: "<<OutputFile<< endl;
  cout<<"setting file: "<<SettingFile<< endl;
  if(HOutputFile == NULL){
    cout<<"No Histogram file created."<<endl;
  }
  int xbin, ybin;

  TFile* File = new TFile(InputFile);
  TIter next(File->GetListOfKeys());
  TKey *key;
  TH2F *h;
  TH1F *ProjHist;
  TH1F *reb4;
  int i=0;
  int j=-1;

  int yrangemin=400;
  int yrangemax=2000;
  double binwidth=1;

  int xrangemin[8];
  int xrangemax[8];
  double gainmin[8];
  double gainmax[8];
  int rangemin[8];
  int rangemax[8];
  int fitmin[8];
  int fitmax[8];
  ifstream input(SettingFile);  
  input.ignore(10000, '\n');

  for(int k=0;k<8;k++){
    input >> xrangemin[k] >> xrangemax[k] >> gainmin[k] >> gainmax[k] >> rangemin[k] >> rangemax[k] >> fitmin[k] >> fitmax[k];
    input.ignore(1000, '\n');
  }

  double gain =0.;
  int gainsteps = 60;
  
  TH2F *Result = new TH2F("Result","Result", yrangemax-yrangemin, yrangemin, yrangemax, gainsteps, gainmin[4], gainmax[4]);

  TH2F *Tilted = new TH2F("Tilted","Tilted", 1000, 0 ,1000, 2000, 0 ,2000);
  double genergy=0;
  double penergy=0;
  double cont=0;
  bool rej = false;
  Double_t oldpar[5] = {0,0,0,100,0};

  ofstream output(OutputFile);
  output << "cont\tmean\tfwhm\tchi\tgain\tname" << endl;
  gStyle->SetOptStat(1111111);
  while ((key=(TKey*)next())){
    i=0;
    oldpar[3]=100;
    if(strcmp(key->GetClassName(),"TH2F") == 0){
      Name = (char*)key->GetName();
      j=lookup(Name);
      //cout << "j " << j << endl;
      //if(strstr(Name,"Forward_gam_pad_0")){
      //if(j>-1){
      if((j>-1) && (j<4)){
	h=(TH2F*) File->Get(Name);
	for(gain=gainmin[j];gain<gainmax[j];gain+=(gainmax[j]-gainmin[j])/gainsteps){
	  ProjHist = new TH1F(Form("%s_proj_%d",Name,i), Form("%s_proj_%.2f",Name,gain), 4000, 0, 4000);
	  for(ybin=yrangemin;ybin<yrangemax;ybin++){
	    genergy=ProjHist->GetBinCenter(ybin);
	    cont=0;
	    for(xbin=xrangemin[j];xbin<xrangemax[j];xbin++){
	      penergy=h->GetXaxis()->GetBinCenter(xbin)-xrangemin[j];
	      //cout << "yrange " << yrangemin <<"-"<<yrangemax<<endl;
	      //cout << "xrange " << xrangemin[j] <<"-"<<xrangemax[j]<<endl;
	      
	      //cout << "ybin " << ybin << " xbin " << xbin << endl;
	      //cout << "g " << genergy << " p " << penergy << " gain " << gain << endl; 
	      //cout << "genergy-penergy*gain " << genergy-penergy*gain << endl;
	      if(genergy-penergy*gain>400 ){
		  cont+=h->GetBinContent(xbin,h->GetYaxis()->FindBin(genergy-penergy*gain));
		Tilted->SetBinContent(xbin,ybin,h->GetBinContent(xbin,h->GetYaxis()->FindBin(genergy-penergy*gain)));
	       }
	    }
	    ProjHist->SetBinContent(ybin,cont);
	    Result->SetBinContent(Result->GetXaxis()->FindBin(ybin), Result->GetYaxis()->FindBin(gain), cont);
	  }
	  
	  Double_t par[5] = {0};
	  reb4 = (TH1F*)ProjHist->Rebin(4);
	  reb4->GetXaxis()->SetRangeUser(rangemin[j],rangemax[j]); 
	  cout << "bin width " << reb4->GetBinWidth(1) << endl;
	  binwidth = reb4->GetBinWidth(1);

	  TF1 *bgline = new TF1("bgline","pol1",rangemin[j],rangemax[j]);
	  reb4->Fit("bgline");

	  par[0] = bgline->GetParameter(0);
	  par[1] = bgline->GetParameter(1);
	  if(rej){
	    cout << "warning: rejecting bg!" << endl;
	    par[0] = 0.;
	    par[1] = 0.;
	  }
	  par[2] = 20;
	  par[3] = 1300;
	  par[4] = 50;
	  //TF1 *fit = new TF1("fit",fpeaksbgpad,rangemin[j],rangemax[j],2+3);
	  TF1 *fit = new TF1("fit",fonegausbg,rangemin[j],rangemax[j],2+3);
	  fit->SetParameters(par);
	  fit->SetParLimits(0,par[0]*0.8,par[0]*1.2);	  
	  fit->SetParLimits(1,par[1]*0.8,par[1]*1.2);	  
	  ////fit->SetParLimits(2,10,40);	  
	  fit->SetParLimits(2,25*binwidth,400*binwidth);	  
	  fit->SetParLimits(3,fitmin[j],fitmax[j]);	  
	  fit->SetParLimits(4,3,1000);	  
	  fit->SetNpx(100);
	  fit->SetLineColor(3);
	  fit->SetLineWidth(1);
	  reb4->Fit("fit","LR");
	  
	  while( fit->GetParameter(3) < oldpar[3] ){
	    cout << "fit: " << fit->GetParameter(3) << " oldpar: " << oldpar[3] << endl;
	    par[2]=oldpar[2];
	    par[3]=oldpar[3];
	    fit->SetParLimits(3,oldpar[3],oldpar[3]+30);	  
	    par[4]=oldpar[4];
	    reb4->Fit("fit","LR");
	  }
	  
	  output << fit->GetParameter(2)/binwidth << "\t" << fit->GetParameter(3) << "\t" << fit->GetParameter(4)*2.35 << "\t" << fit->GetChisquare() << "\t" << gain << "\t" << Name <<"_"<< i << endl;
	  TFile* HFile = new TFile(HOutputFile,"update");
	  reb4->Write("",TObject::kOverwrite);
	  fit->Write("",TObject::kOverwrite);
	  bgline->Write("",TObject::kOverwrite);
	  Tilted->Write(Form("%s_titled_%d",Name,i),TObject::kOverwrite);
	  Result->Write("",TObject::kOverwrite);
	  HFile->Close();
	  i++;
	  for(int k=0;k<5;k++)
	    oldpar[k]=fit->GetParameter(k);
	}
      }
    }
  } 
  return 0;
}
